The dissociation energy of N2 
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Abstract . - . - . 

We examine the requirements for very accurate ab initio quantum chemical 
prediction of dissociation energies, using a detailed investigation of the nitrogen 
molecule. Although agreement with experiment to within 1 kcal/mol is not achieved 
even with the most elaborate multireference CI fl wave functions and largest basis sets 
currently feasible, it is possible to obtain agreement to within about 2 kcal/mol, or 
1% of the dissociation energy. At this level it is necessary to account for core- valence 
correlation effects and to include up to /i-type functions in the basis. We have also 
investigated the effect of f-type functions, the use of different reference configuration 
spaces, and basis set superposition error,. After discussing these results, we speculate 
on the remaining sources of error in our best calculations^ 


f Mailing address: NASA Ames Research Center, Moffett Field, CA 94035 
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I. Introduction 

The dissociation energy ( D e ) of the nitrogen molecule has often been used as 
a test of the capabilities of quantum chemical methods. The experimental value 
of 228.4 kcal/mol [1] is larger than most diatomic D e values, and there is a very 
large correlation contribution, with the SCF result being little more than half that 
of experiment. An accurate computed D e value thus requires a very sophisticated 
treatment of the correlation problem. (Our discussion here excludes calculations 
in which good agreement with experiment derives from cancellations of errors or 
from basis set superposition error). In a recent study using an elaborate multi- 
reference Cl (MRCI) wave function and a basis set with up to #-type functions, 
LanghofF, Bauschlicher and Taylor (LBT) [2] obtained a D e value of 224 kcal/mol. 
This is the same as the value obtained by Binkley and Frisch [3] using an MP4 
treatment of correlation and a smaller basis set: the comparison of single-reference 
results given by Ahlrichs and co-workers [4] (who obtained 214 kcal/mol using the 
coupled-pair functional approach and a basis set somewhat smaller than Ref. 2) 
suggests that MP4 probably overestimates D e . However, the LBT result is itself 
actually an overestimate, by some 2.5 kcal/mol, because the deletion of cr virtual 
orbitals necessitated by near linear dependence degrades the asymptotic description 
by this amount. Further, correction for basis set superposition error would reduce 
the LBT result by almost 0.5 kcal/mol, and the D e value of Binkley and Frisch by 
almost 3 kcal/mol. Hence the best ab initio estimate so far, in the sense of degree of 
correlation treatment and size and saturation of basis set, is around 221 kcal/mol. 

An error of 7.5 kcal/mol in a D e value of 228.4 kcal/mol is not a large relative 
error, but it is considerably larger than the commonly-stated target of “chemical 
accuracy” (1 kcal/mol). If the inadequacies of quantum chemical methodology 
are to be understood, it would be desirable to investigate how the N 2 ground- 
state D e value can be computed to higher accuracy. A number of possible sources 
of error can be considered. First, the correlation treatment could be inadequate. 
This might involve either inadequacies in the MRCI reference space, or the lack 
of size-consistency in the Cl expansion. Further, the LBT results were obtained 
from a treatment in which the core electrons were not correlated, and the core 
correlation results of Ahlrichs et al [4] suggest that this could account for almost 
1 kcal/mol. Second, the basis set used by LBT includes up to g-type functions: 
there may thus be errors from neglect of higher angular functions. There may also 
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be errors from lack of saturation of the primitive or contracted sets used. Third, 
the question of basis set superposition error, which may give a spurious increase in 
D e , should be investigated. Finally, we may note here that relativistic effects, as 
estimated by first-order perturbation theory [5], contribute less than 0.1 kcal/mol 
to the binding (actually decreasing D e very slightly). They can be neglected for our 
present purposes. Spin-orbit effects should be negligible for the N 2 ground-state 
dissociation as all the atomic 2 p orbitals participate in the bonding. It also seems 
reasonable to assume that other effects such as Born-Oppenheimer breakdown terms 
can also be neglected. 

In the present work we shall discuss in detail all of the above potential sources 
of error. A brief description of the computational methods is given in the next 
section, but it will be more convenient to discuss the details of wave functions and 
basis sets in the later sections. In Section III we consider the effects of extending 
the correlation treatment, including the effects of core-valence correlation and the 
question of size-consistency. The effect of improving the basis set, including both 
extending the primitive and contracted s, p, d , /, and g spaces and adding h- and 
i-type functions, is described in Section IV, together with the question of basis 
set superposition error. Our best theoretical estimates of D e and other spectro- 
scopic constants, together with possible sources of remaining errors, are discussed 
in Section V, and our conclusions are given in Section VI. 

II. Computational methods. 

The basis sets used in the present calculations are derived from two dif- 
ferent primitive Gaussian basis sets: a (13s 8p 6d 4/ 2 g) set and a larger 

(18s 13p 6d 5/ 4 <7 3 h 2 i) set. The latter should be close to saturation in all / values 
included. The sp sets are from van Duijneveldt [6] and Partridge [7], respectively. 
The polarization sets are even-tempered with a ratio of 2.5 between successive ex- 
ponents; these are based on geometric mean exponents of (1.2) ,_2 . The primitive 
sets were contracted using atomic natural orbitals [8] for the N ( 4 5) ground state. 
The various contracted sets used are described in Sections III and IV . 

The molecular orbitals used were obtained from CASSCF wave functions in 
which the 2p- derived orbitals and electrons were active, giving six active electrons 
in six orbitals in the CASSCF calculations. The use of orbitals from a larger active 
space CASSCF calculation, in which all ten valence electrons are active in the 
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eight valence orbitals, was investigated, but near r e the two active spaces yield 
almost the same orbitals, while at laxge internucleax distances the larger active 
space creates difficulties in defining the orbitals for a subsequent Cl calculation 
because of Is /2s orbital mixing. We therefore report here only results obtained 
with molecular orbitals from the smaller CASSCF active space. 

All MRCI wave functions used in this work are based on CAS reference spaces, 
that is, all configuration state functions (CSFs) that can be generated with a par- 
ticular choice of active space are used as references. Two active spaces have been 
used, one in which six electrons are distributed among six active orbitals, nominally 
derived from the atom 2 p orbitals, and a larger space in which all ten valence elec- 
trons are distributed among the eight valence orbitals. The calculations are labelled 
MRCI(6act) and MRCI(lOact), respectively. It should be noted that ten electrons 
are correlated in all MRCI wave functions, the difference arises from whether the 
2s electrons and orbitals are included in the active space or not. In addition to the 
MRCI results, we also quote results obtained using a multireference analog [9] of 
Davidson’s correction [10], and with the averaged coupled-pair functional (ACPF) 
method of Gdanitz and Ahlrichs [11]. 

Spectroscopic constants (r e , u e , and D e ) were obtained by fitting computed 
energy points to a polynomial in 1 /r, as described in detail in Section V below. 
The use of 1/r rather than r as the independent variable provides a more rapidly 
convergent expansion of the potential energy about a given point [12]. For the 
investigation of changes in basis sets and correlation treatment the “dissociation 
energy” was computed as the difference between the total energy at 2.1 ao (close to 
the experimental equilibrium) and 100 «o • Tests show that the differential effects 
associated with a change in the level of treatment would not be affected by the use 
of another (similar) r value. 

All calculations were performed using the MOLECULE-SWEDEN [13,14] pro- 
gram system on the CRAY X— MP/48 and NAS Facility CRAY Y— MP/832 at 
NASA Ames and the CRAY-2 at the Minnesota Supercomputer Center. 

III. Convergence of D e with correlation treatment. 

The D e values obtained from different correlation treatments in a 
[5s 4 p 3d 2/ 1 g] basis (contracted from the smaller primitive set) are listed in 
Table 1. The MRCI(6act) calculations correspond to those of LBT, although as 
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noted in the introduction their result was somewhat larger because of an unbalanced 
treatment of the atoms relative to the molecule. The close agreement between the 
MRCI(6act)+Q and ACPF(6act) results suggest that the multireference Davidson 
correction for higher excitations gives reliable results here. The D e value obtained 
from the MRCI(lOact) calculation is very similar to that from MRCI(6act)+Q (or 
ACPF(6act)), but the effect of the +Q correction (or of ACPF) on the MRCI(lOact) 
calculation is to decrease D e , which is somewhat unexpected. The different behav- 
ior with the different MRCI reference spaces can be understood by examining their 
performance on the states of nitrogen atom. 

Table 2 lists total energies for the lowest states of nitrogen atom (derived from 
the occupation 2s 2 2p 3 ) and separations between them. The separation between the 
two lowest states, 4 5 and 2 £>, is not affected by the change from an SCF-based 
description (singles and doubles Cl (SDCI) equivalent to a three active orbital 
CASSCF reference or to the 6act treatment of N 2 ) to an MRCI based on a five- 
electron four active orbital CASSCF treatment (equivalent to the lOact treatment of 
N 2 ), because only one CSF of the appropriate symmetry can be formed in this active 
space in the atom. This is not the case for the i S — 2 P separation, as the CASSCF 
treatment for the atom accounts for 2s /2p near- degeneracy mixing in the 2 P state. 
As the 4 S'- 2 P separation is overestimated by the 6act treatment of N 2 , the mixing 
of configurations derived from the 2 P + 2 P limit into the molecular wave function 
will be underestimated, reducing the strength of the bond and underestimating D e . 
Hence it is reasonable that the MRCI(lOact) treatment produces a larger D e than 
MRCI(6act). Incidentally, the discrepancies between the best calculated atomic 
term separations given in Table 2 and the experimental values [15] are entirely due 
to basis set deficiencies: the use of a [6s 5 p 4 d 3/ 2 g 1 h li] ANO basis reduces 
the separations by about 0.04 eV at the correlated level, bringing them into almost 
perfect agreement with experiment. 

The rather unexpected reduction in D e on going from the MRCI(lOact) treat- 
ment to MRCI(10act)+Q can also be understood in terms of N atom. The most 
important correlating configuration in the 4 S ground state is the single orbital exci- 
tation 2s — * 3d, accompanied by a recoupling of the 2 p electrons. This is essentially 
a purely atomic effect, as we have noted before [2]. For the MRCI(6act) calcula- 
tion at infinite separation this excitation can occur on only one of the atoms at a 
time, because of the spin and symmetry coupling constraints on the reference CSFs. 
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However, in the MRCI(lOact) calculation this effect can occur on both atoms at 
the same time because there are far fewer coupling restrictions on the active oc- 
cupations compared to the MRCI(6act) case. The possibility of recovering more 
atomic correlation at the MRCI(lOact) level results in a reduction in the weight of 
the reference CSFs in the MRCI(lOact) wave function relative to MRCI(6act), at 
infinite separation, while near r e the larger active space calculation has an increased 
weight of the reference CSFs relative to MRCI(6act). Of course, as the larger cal- 
culation allows important atomic correlation effects on both atoms simultaneously 
it can legitimately be regarded as more size-consistent than MRCI(6act), and so 
inclusion of the +Q correction results in double counting some effects. We conclude 
that while the MRCI(6act)+Q results are to be preferred over MRCI(6act), the 
uncorrected MRCI(lOact) results should be more reliable than MRCI(10act)+Q. 

An effect excluded from the calculations discussed so far is core correlation. In 
CPF calculations in which all 14 electrons were correlated, Ahlrichs and co-workers 
obtained a increase of some 0.7 kcal/mol in D e as a result of core correlation [4]. 
Experience suggests that the core-core correlation will be unimportant and that 
this effect represents core- valence correlation. We have computed the effect of core- 
valence correlation on D e by performing multireference Cl calculations based on 
MRCI(6act), but with additional CSFs in which one electron is excited from the 
core orbitals together with one (or none) from the valence orbitals. The basis set 
used for these calculations comprised the uncontracted (13s 8 p 6d) set described 
above, augmented with the / ANO with the largest occupation number from the 
smaller ANO set. Tests on core- valence correlation in nitrogen atom demonstrated 
that higher-exponent d functions or other higher angular momentum functions made 
a negligible contribution. The core-valence contribution to D e was computed to be 
0.7 kcal/mol, in remarkable agreement with the work of Ahlrichs et al [4]. From the 
lack of sensitivity of the atomic core- valence correlation to changes in basis set or 
correlation treatment it can safely be concluded that the D e value of N 2 is increased 
by less than 1 kcal/mol as a result of core- valence correlation. 

IV. Convergence of D e with basis set. 

Table 3 shows the variation in the D e value of N 2 with atomic basis set. 
These results were obtained at the MRCI(6act) level. (The difference between 
MRCI(lOact) and MRCI(6act) was explored for several different basis sets, but was 
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found to be the same (0.7 kcal/mol in D e ) for all the large basis sets.) We consider 
first the effect of expanding the contracted set obtained from the (13s 8 p 6 d 4/ 2 g) 
primitive set. The extension from [4s 3 p 2d If] to [5s 4 p 3 d 2 /] and then to 
[6s 5 p 4 d 3/] indicates that the contraction error in D e in the latter basis should be 
less than 0.5 kcal/mol. Adding a diffuse sp set increases D e by only 0.2 kcal/mol. 
Comparing the [5s 4p 3d 2/ 1 g] and [5s 4 p 3 d 2 /] results shows that the g ANO 
increases D e by 2.5 kcal/mol: this is a little smaller than the 3.0 kcal/mol observed 
by Ahlrichs et al [4], but the smaller primitive sets used in Ref. 4 may lead to ef- 
fects which are actually due to d and / set unsaturation being attributed to the g 
set. Adding spdfg ANOs to the [5s 4 p 3 d 2/ 1 g] set increases D e by 1.5 kcal/mol: 
comparing this number with the 1.2 kcal/mol observed on adding spdf ANOs to 
the [5s 4p 3d 2/] set suggests that additional g orbitals would increase D e by only 
0.3 kcal/mol. Finally, adding a single h primitive set to the [6s 5p 4d 3/ 2 g] basis 
increases D e by 0.7 kcal/mol, suggesting that the total contribution of h and higher 
angular momentum functions will be on the borderline of “chemical accuracy” . 

The effect of replacing the (13s 8p 6d 4/ 2 g) primitive set with the 
(18s 13p 6d 5/ 4 g) primitive set is to increase D e by only 0.1 kcal/mol in a 
[5s 4p 3d 2/ 1 g] contracted set. Thus the smaller primitive set seems to be sat- 
urated in all the angular momentum quantum numbers, at least as far as the D e 
value is concerned. Adding an h - type ANO, contracted from 3 primitive functions, 
increases D e by 1.1 kcal/mol, but this is actually an overestimate of the h effect, as 
there is partial compensation for unsaturation of the spdfg space. When the h ANO 
is added to the [6s 5 p 4d 3/ 2 g] ANO set the increase in D c is 0.8 kcal/mol, only 
slightly more than the effect of the single h primitive, and consistent with the effect 
of increasing the contracted set from [5s 4p 3d 2/ 1 g 1 h] to [6s 5p 4d 3/ 2 g 1 h]. 
This gives an increase in D e of 1.1 kcal/mol. Finally, the addition of an i-type 
ANO, contracted from two primitive i functions, increases D e by 0.4 kcal/mol. It 
therefore seems safe to conclude that the contribution to D e from higher than i-type 
functions will be less than 0.3 kcal/mol. 

The need to include g- and h - type functions in the basis set to achieve chem- 
ical accuracy in the N 2 value creates rather large basis sets. It is useful to 
examine the contributions of the different components of the angular sets, in case 
some are noticeably more important than others. If this is the case, the basis set 
could be reduced by eliminating the unimportant components, with a consequent 
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reduction in the overall computational effort. This possibility has been investigated 
by computing the contribution to D e made by the /, g and h components in the 
[ 5 s 4 p 3 d 2 / 1 g l/i] basis (using the small primitive set for the / and g inves- 
tigations and the large set for the h case). Various approaches were tried, such 
as selectively deleting each angular component or adding only the desired angular 
component, but all gave essentially the same results. Table 4 displays the increase 
in D e obtained with successive addition of different components. First, it appears 
that there is little scope for reducing the basis set, as the differences between the 
contribution from different angular components is seldom large. Further, the main 
contribution in the / case is to the tv and 6 spaces, while in the g case this shifts to 
the 8 and <f> spaces. This is perhaps surprising: it seems plausible to expect the 
functions to be the most important when g functions are first added, for instance, 
as they are the first functions of 7 type to be included. However, this is not what is 
observed. It is also clear the the contribution of higher angular functions becomes 
more “isotropic” as the angular quantum number increases, as the h example in 
Table 4 shows. This effect becomes more pronounced on noting that in the case 
of tv or higher components two members of a degenerate pair of orbitals are added 
together. Overall, then, it must be concluded that it is necessary to include all 
members of any set of higher angular functions that are added. 

The last aspect of our basis set investigations is the question of basis set su- 
perposition error (BSSE). In the smaller primitive basis the BSSE contribution to 
D e (that is, twice the energy lowering obtained by including a ghost basis in an 
atom calculation) is computed to be 0.28 kcal/mol in the [ 6 s 5 p 4 d 3/ 2 g l/i] con- 
tracted set, to be compared with 0.36 and 0.73 kcal/mol in the [ 6 s 5 p 4 d 3/ 2 g] and 
[ 6 s 5 p 4 d 3/] sets, respectively. This is determined using a counterpoise correction 
with the full ghost basis. Thus in our largest calculations BSSE should not exceed 
0.3 kcal/mol. 

V. Computed spectroscopic constants. 

Table 5 contains our best directly computed results for D e , u; e , and r e . These 
results are obtained from a sixth-order fit in 1 /r to eleven points around equilibrium, 
obtained from MRCI( 6 act) calculations in a [ 6 s 5p 4d 3/ 2 g 1 h] contracted basis 
(large primitive set). The agreement with the experimentally derived values for r e 
and u e is very good, but the computed D e value is still 3.2 kcal/mol less than the 
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experimental value. From the results presented in the preceding sections we can 
add 0.7 kcal/mol to our computed value to account for core-valence correlation, add 
0.4 kcal/mol for the effect of adding an f-type function, and subtract 0.3 kcal/mol 
for the effects of BSSE. Thus based on the various calculations presented above we 
obtain an estimate of 226.0 kcal/mol, which is still outside chemical accuracy in 
terms of agreement with experiment, but is accurate to essentially 1% of D e . 

We may speculate on the remaining sources of error in the theoretical D e value. 
By extrapolating from the basis set investigations of Section IV it seems very likely 
that additional spdfghi functions would add at most 0.5 kcal/mol to D e , and the 
contribution of higher than 7 -type functions can hardly exceed 0.3 kcal/mol. These 
effects would reduce the discrepancy to 1.6 kcal/mol. It seems unlikely that addi- 
tional core or core- valence correlation could give rise to an effect of this magnitude, 
and it is highly improbable that more exotic phenomena, such as Born-Oppenheimer 
breakdown or the diagonal Born-Oppenheimer correction, could contribute to even 
within an order of magnitude. Relativistic effects have also been shown to be neg- 
ligible, as discussed above. It therefore seems likely that the valence correlation 
treatment is still in error by 1.6 kcal/mol. This is not inconsistent with the molec- 
ular consequences of quite subtle changes in the description of the atomic limits, as 
discussed in Section III above. The best way to resolve this question would be to 
perform 10-electron full Cl calculations on N 2 , using a realistic basis set such as the 
valence double zeta plus polarization set used by Bauschlicher and Langhoff [16]. 
Recent advances in the computational methodology of full Cl calculations [17,18] 
should make such a calculation possible in the near future. 

VI. Conclusions. 

Even with very elaborate MRCI wave functions and extended ANO basis sets 
with high angular momentum functions, it is not yet possible to compute the 
dissociation energy of N 2 from first priniciples to within 1 kcal/mol — a some- 
what arbitrary definition of chemical accuracy. However, an accuracy of about 1% 
does appear to be achievable. From the results of this work it appears safe to as- 
sert that the basis set limit for a calculation using the MRCI(6act) configuration 
space and including Davidson’s correction for higher excitations is a D e value of 
226.8 kcal/mol. The remaining 1.6 kcal/mol difference with experiment is most 
probably attributable to inadequacies in the valence correlation treatment, suggest- 
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ing that larger reference spaces will be required. 
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Table 1. Convergence of N 2 D e (kcal/mol) with correlation treatment 0 . 


Method 

6 active electrons 

10 active electrons 

Reference 6 

204.5 

211.9 

MRCI 

222.1 

222.8 

MRCI+Q 

222.6 

221.6 

ACPF 

222.5 

221.5 

Ref. weight 0 (2.1 do) 

93.7 

94.2 

Ref. weight (100 do) 

93.3 

93.0 


a [55 4 p 3d 2/ lg] ANO basis. 

6 Energy obtained with reference CSFs. 

c Weight of reference CSFs in MRCI wave function (percentage) at stated r(NN). 
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Table 2. N atom term separations (eV)“. 


Method 

4 S- 2 D 

4 S — 2 P 

SCF 

2.85 

4.71 

SDCI 

2.45 

3.72 

SDCI+Q 

2.42 

3.58 

CASSCF 6 


3.77 

MRCI 


3.65 

MRCI+Q 


3.64 

Expt c 

2.38 

3.58 


° [5s 4 p 3d 2/ lg] ANO basis. 

b 5 active electrons in 4 active orbitals. This reduces to an SCF wave function for 
the 4 S' and 2 D states. 
c Moore, Ref. 14. 


13 



Table 3. Convergence of D e (kcal/mol) with atomic basis set®. 


Basis set 


D e 


(13s 8p 6d 4/ 2 g 1 h) primitive set 


[4s 3 p 2d 1/] 

216.0 

[5s 4p 3d 2/] 

219.6 

[6s 5 p 4 d 3/] 

220.8 

[6s 5 p 4 d 3/]+ diffuse sp 

221.0 

[5s 4p 3d 2/ 1 g\ 

222.1 

[6s 5p 4d 3/ 2 g] 

223.6 

[6s 5p 4d 3/ 2 g l/i] 

224.3 

(18s 

13p 6 d 5/ 4 g 3 h 2 i) primitive set 

[5s 4p 3d 2/ 1 g\ 

222.3 

[5s 4p 3d 2/ 1 g 1 h] 

223.4 

[5s 4p 3d 2/ 1 g 1 h If] 

224.0 

[6s 5 p 4d 3/ 2 g] 

223.7 

[6s 5 p 4d 3/ 2<? l/i] 

224.5 

[6s 5 p 4 d 3} 2g lh If] 

225.0 


a MRCI(6act) calculation 
6 See text for details of primitive basis sets. 
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Table 4. D e contributions (kcal/mol) from different angular components®. 


Component 

2/ set 

1 g set 

lh set 

a 

0.81 

0.19 

0.10 

7 r 

2.98 

0.57 

0.11 

6 

3.08 

0.72 

0.16 


1.39 

0.71 

0.20 

7 


0.32 

0.15 

7 



0.06 


° MRCI(6act) calculation; small primitive basis for / and g case, added to a [5s 4 p 3d] 
ANO set. Large primitive basis for h case, added to a [6s 5p 4 d] ANO set. 
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Table 5. Spectroscopic constants for N2 a . 



r e (A) 

uj e (cm 1 ) 

£> e (kcal/mol) 

Computed 

1.100 

2353 

225.2 

Experiment 6 

1.098 

2359 

228.4 


a MRCI(6act)+Q calculation; [65 5 p 4 d 3/ 2 g lft] ANO basis (large primitive set). 
6 Huber and Herzberg, Ref. 1. 
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